{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4ed546c8-cbad-4721-ba6b-7d185222990e",
   "metadata": {},
   "source": [
    "# Example 4\n",
    "No source, Multiply-connected domain, Dirichlet/Neumann problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "228a690d-545c-48b5-a2d5-6a3e40d7963b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import calfem.geometry as cfg\n",
    "import calfem.mesh as cfm\n",
    "import calfem.vis as cfv\n",
    "import numpy as np\n",
    "from toolkits import gauss_line_int, gauss_tri_int\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b47a234-100f-4115-a5ee-670e2b590fdd",
   "metadata": {},
   "source": [
    "## 1. Geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1002c549-ae11-4ae3-b712-f909f8ba093c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Number of holes\n",
    "nH = 1\n",
    "# Mesh size factor\n",
    "factor = 20\n",
    "\n",
    "# externel boundary ID\n",
    "dir_id_1 = 100\n",
    "dir_id_2 = 200\n",
    "dir_id_3 = 300\n",
    "dir_id_4 = 400\n",
    "# hole boundaryd ID\n",
    "hole_neu_id_1 = 500\n",
    "hole_neu_id_2 = 600\n",
    "hole_neu_id_3 = 700\n",
    "hole_neu_id_4 = 800\n",
    "\n",
    "# Geometry definition\n",
    "g = cfg.Geometry()\n",
    "g.point([0.0, 0.0], ID=0)\n",
    "g.point([1.0, 0.0], ID=1)\n",
    "g.point([1.0, 1.0], ID=2)\n",
    "g.point([0.0, 1.0], ID=3)\n",
    "g.point([0.4, 0.4], ID=4)\n",
    "g.point([0.6, 0.4], ID=5)\n",
    "g.point([0.6, 0.6], ID=6)\n",
    "g.point([0.4, 0.6], ID=7)\n",
    "g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "g.spline([1, 2], el_on_curve=factor, marker=dir_id_2, ID=1) \n",
    "g.spline([2, 3], el_on_curve=factor, marker=dir_id_3, ID=2) \n",
    "g.spline([3, 0], el_on_curve=factor, marker=dir_id_4, ID=3)\n",
    "g.spline([4, 5], el_on_curve=factor/2, marker=hole_neu_id_1, ID=4)\n",
    "g.spline([5, 6], el_on_curve=factor/2, marker=hole_neu_id_2, ID=5)\n",
    "g.spline([6, 7], el_on_curve=factor/2, marker=hole_neu_id_3, ID=6)\n",
    "g.spline([7, 4], el_on_curve=factor/2, marker=hole_neu_id_4, ID=7)\n",
    "g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "cfv.drawGeometry(g)\n",
    "cfv.showAndWait()\n",
    "\n",
    "# g.point([0.0, 0.0], ID=0)\n",
    "# g.point([1.0, 0.0], ID=1)\n",
    "# g.point([1.0, 1.0], ID=2)\n",
    "# g.point([0.0, 1.0], ID=3)\n",
    "# g.point([0.5, 0.5], ID=4)\n",
    "# g.point([0.6, 0.5], ID=5)\n",
    "# g.point([0.5, 0.6], ID=6)\n",
    "# g.point([0.4, 0.5], ID=7)\n",
    "# g.point([0.5, 0.4], ID=8)\n",
    "# g.spline([0, 1], el_on_curve=factor, marker=dir_id_1, ID=0)\n",
    "# g.spline([1, 2], el_on_curve=factor, marker=neu_id_1, ID=1) \n",
    "# g.spline([2, 3], el_on_curve=factor, marker=dir_id_2, ID=2) \n",
    "# g.spline([3, 0], el_on_curve=factor, marker=dir_id_3, ID=3)\n",
    "# g.circle([5, 4, 6], el_on_curve=factor/4, ID=4)\n",
    "# g.circle([6, 4, 7], el_on_curve=factor/4, ID=5)\n",
    "# g.circle([7, 4, 8], el_on_curve=factor/4, ID=6)\n",
    "# g.circle([8, 4, 5], el_on_curve=factor/4, ID=7)\n",
    "# g.surface([0, 1, 2, 3],[[4,5,6,7]])\n",
    "# # g.structuredSurface([0, 1, 2, 3])\n",
    "# Display mesh\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "90c4263d-6fe7-484d-8ca4-8a590e25a3b6",
   "metadata": {},
   "source": [
    "## 2. Meshing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "07931634-1ddc-4c54-9b23-e70b883c5d48",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Info    : GMSH -> Python-module\n",
      "[[  4 119]\n",
      " [119 118]\n",
      " [118 117]\n",
      " [117 116]\n",
      " [116 115]\n",
      " [115 114]\n",
      " [114 113]\n",
      " [113 112]\n",
      " [112 111]\n",
      " [111   7]\n",
      " [  7 110]\n",
      " [110 109]\n",
      " [109 108]\n",
      " [108 107]\n",
      " [107 106]\n",
      " [106 105]\n",
      " [105 104]\n",
      " [104 103]\n",
      " [103 102]\n",
      " [102   6]\n",
      " [  6 101]\n",
      " [101 100]\n",
      " [100  99]\n",
      " [ 99  98]\n",
      " [ 98  97]\n",
      " [ 97  96]\n",
      " [ 96  95]\n",
      " [ 95  94]\n",
      " [ 94  93]\n",
      " [ 93   5]\n",
      " [  5  92]\n",
      " [ 92  91]\n",
      " [ 91  90]\n",
      " [ 90  89]\n",
      " [ 89  88]\n",
      " [ 88  87]\n",
      " [ 87  86]\n",
      " [ 86  85]\n",
      " [ 85  84]\n",
      " [ 84   4]]\n",
      "[[ 5 92]\n",
      " [92 91]\n",
      " [91 90]\n",
      " [90 89]\n",
      " [89 88]\n",
      " [88 87]\n",
      " [87 86]\n",
      " [86 85]\n",
      " [85 84]\n",
      " [84  4]]\n",
      "[[  6 101]\n",
      " [101 100]\n",
      " [100  99]\n",
      " [ 99  98]\n",
      " [ 98  97]\n",
      " [ 97  96]\n",
      " [ 96  95]\n",
      " [ 95  94]\n",
      " [ 94  93]\n",
      " [ 93   5]]\n",
      "[[  7 110]\n",
      " [110 109]\n",
      " [109 108]\n",
      " [108 107]\n",
      " [107 106]\n",
      " [106 105]\n",
      " [105 104]\n",
      " [104 103]\n",
      " [103 102]\n",
      " [102   6]]\n",
      "[[  4 119]\n",
      " [119 118]\n",
      " [118 117]\n",
      " [117 116]\n",
      " [116 115]\n",
      " [115 114]\n",
      " [114 113]\n",
      " [113 112]\n",
      " [112 111]\n",
      " [111   7]]\n",
      "Dirichlet !\n"
     ]
    }
   ],
   "source": [
    "mesh = cfm.GmshMesh(g,return_boundary_elements=True)\n",
    "mesh.elType = 2          # Degrees of freedom per node.\n",
    "mesh.dofsPerNode = 1     # Factor that changes element sizes.\n",
    "# mesh.elSizeFactor = 0.5 # Element size Factor\n",
    "\n",
    "coords, edof, dofs, bdofs, elementmarkers, boundaryElements  = mesh.create()\n",
    "\n",
    "#############################\n",
    "#                           #\n",
    "# Elements, Edges and Nodes #\n",
    "#                           #\n",
    "#############################\n",
    "element = edof - 1\n",
    "node = coords\n",
    "edge = dict()\n",
    "for i, e in enumerate(element):\n",
    "    edge[(e[0], e[1])] = i\n",
    "    edge[(e[1], e[2])] = i\n",
    "    edge[(e[2], e[0])] = i\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#   Externel edges   #\n",
    "#                    #\n",
    "######################\n",
    "dir_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_1]])-1\n",
    "dir_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_2]])-1\n",
    "dir_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_3]])-1\n",
    "dir_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[dir_id_4]])-1\n",
    "\n",
    "dir_ids = [dir_id_1, dir_id_2, dir_id_3, dir_id_4]\n",
    "dir_edges = set([tuple(np.array(edge['node-number-list'])-1) for marker in dir_ids for edge in boundaryElements[marker]])\n",
    "neu_edges = set([])\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#     Hole edges     #\n",
    "#                    #\n",
    "######################\n",
    "# first hole\n",
    "hole_neu_edge_1 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_1]])-1\n",
    "hole_neu_edge_2 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_2]])-1\n",
    "hole_neu_edge_3 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_3]])-1\n",
    "hole_neu_edge_4 = np.array([edge['node-number-list'] for edge in boundaryElements[hole_neu_id_4]])-1\n",
    "\n",
    "hole_neu_edges = [hole_neu_edge_1, hole_neu_edge_2, hole_neu_edge_3, hole_neu_edge_4]\n",
    "\n",
    "ccw_hole_neu_edges = np.array([edge[::-1] for edges in hole_neu_edges[::-1] for edge in edges[::-1]])\n",
    "\n",
    "ccw_hole_neu_edge_1 = np.array([edge[::-1] for edge in hole_neu_edge_1[::-1]])\n",
    "ccw_hole_neu_edge_2 = np.array([edge[::-1] for edge in hole_neu_edge_2[::-1]])\n",
    "ccw_hole_neu_edge_3 = np.array([edge[::-1] for edge in hole_neu_edge_3[::-1]])\n",
    "ccw_hole_neu_edge_4 = np.array([edge[::-1] for edge in hole_neu_edge_4[::-1]])\n",
    "\n",
    "print(ccw_hole_neu_edges)\n",
    "print(ccw_hole_neu_edge_1)\n",
    "print(ccw_hole_neu_edge_2)\n",
    "print(ccw_hole_neu_edge_3)\n",
    "print(ccw_hole_neu_edge_4)\n",
    "\n",
    "#################################\n",
    "#                               #\n",
    "#  Crossing elements and edges  #\n",
    "#                               #\n",
    "#################################\n",
    "direction = np.array([0, -1])\n",
    "Bh = []\n",
    "Eh = {}\n",
    "\n",
    "Bh.append(tuple(hole_neu_edge_1[0]))\n",
    "\n",
    "# crossing flag \n",
    "connected = False\n",
    "\n",
    "######################################################\n",
    "# TODO: New structures of Eh, Bh for multiple holes. #\n",
    "######################################################\n",
    "while not connected:\n",
    "    cross_edge = Bh[-1]\n",
    "    # print(cross_edge)\n",
    "    # print(node[cross_edge[0]])\n",
    "    # print(node[cross_edge[1]])\n",
    "    if (cross_edge not in dir_edges) and (cross_edge not in neu_edges):\n",
    "        cross_edge = cross_edge[::-1]\n",
    "        ei = edge[cross_edge]\n",
    "        \n",
    "        i, j, k = element[ei]\n",
    "        order = [i,j,k]  \n",
    "        \n",
    "        local_edges = set([(i,j), (j,k), (k,i)])\n",
    "        local_edges.remove(cross_edge)\n",
    "        e1 = np.array(local_edges.pop())\n",
    "        e2 = np.array(local_edges.pop())\n",
    "        \n",
    "        e1n = node[e1[1]] - node[e1[0]]\n",
    "        e1n = e1n/np.linalg.norm(e1n)\n",
    "        \n",
    "        e2n = node[e2[1]] - node[e2[0]]\n",
    "        e2n = e2n/np.linalg.norm(e2n)\n",
    "        \n",
    "        d1 = np.abs(np.dot(direction, e1n))\n",
    "        d2 = np.abs(np.dot(direction, e2n))\n",
    "        if d1 < d2:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e1[0]] - node[e1[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e1[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e1))\n",
    "        else:\n",
    "            ae = np.zeros(3)\n",
    "            L1 = np.linalg.norm(node[cross_edge[0]] - node[cross_edge[1]])\n",
    "            L2 = np.linalg.norm(node[e2[0]] - node[e2[1]])\n",
    "            ae[order.index(cross_edge[0])]= 1/L1\n",
    "            ae[order.index(e2[0])]= -1/L2\n",
    "            Eh[ei] = ae\n",
    "            Bh.append(tuple(e2))\n",
    "            \n",
    "        # print(e1,'\\t',e1n)\n",
    "        # print(e2,'\\t',e2n)\n",
    "        # print(ae)\n",
    "        # print(d1,d2)\n",
    "        # print(\"----------------------------------\")\n",
    "        \n",
    "            \n",
    "    else:\n",
    "        if cross_edge in neu_edges:\n",
    "            print('Neumann ! ')\n",
    "        else:\n",
    "            print('Dirichlet !')\n",
    "        connected = True\n",
    "    \n",
    "# print(Eh)\n",
    "        \n",
    "# np.all()\n",
    "\n",
    "######################\n",
    "#                    #\n",
    "#  Numbers of sets   #\n",
    "#                    #\n",
    "######################\n",
    "E = len(element)\n",
    "N = len(node)\n",
    "num_dof = N + nH\n",
    "\n",
    "# Display mesh\n",
    "cfv.draw_mesh(coords=coords, edof=edof, dofs_per_node=mesh.dofs_per_node,\n",
    "              el_type=mesh.el_type,filled=True,title=\"Example 2 - Mesh\")\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edf91826-efaf-416e-9306-d783aa7ea2c6",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Boundary condition"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cb26e870-b29d-41b5-bb21-0c83e022121b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Externel boundary condition\n",
    "########## Dirichlet ##########\n",
    "def dir1(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir2(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir3(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "def dir4(x, y):\n",
    "    return np.zeros_like(x)\n",
    "########## Neumann ############\n",
    "# ...\n",
    "\n",
    "# Hole's boundary condition\n",
    "########## Dirichlet ##########\n",
    "# ...\n",
    "########## Neumann ############\n",
    "def hole_neu1(x, y):\n",
    "    return (np.sin(x)/np.sin(1))*(np.cosh(0.4)/np.sinh(1))\n",
    "\n",
    "def hole_neu2(x, y):\n",
    "    return -(np.cos(0.6)/np.sin(1))*(np.sinh(y)/np.sinh(1))\n",
    "\n",
    "def hole_neu3(x, y):\n",
    "    return -(np.sin(x)/np.sin(1))*(np.cosh(0.6)/np.sinh(1))\n",
    "\n",
    "def hole_neu4(x, y):\n",
    "    return (np.cos(0.4)/np.sin(1))*(np.sinh(y)/np.sinh(1))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d5fe1912-f8d3-4895-9e96-3156355eb71c",
   "metadata": {},
   "source": [
    "## 3. Assemble flexibility matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "fb18dbeb-89fc-4566-b9bc-7627f4bef33f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "size: 922 x 922\n",
      "rank: 921\n"
     ]
    }
   ],
   "source": [
    "F = np.zeros((num_dof, num_dof))\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    # print(e)\n",
    "    # print(node.shape)\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]])\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        # print(Eh[i])\n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        ### Todo: multiple holes ###\n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        F[np.ix_(ey,ey)] = F[np.ix_(ey,ey)] + F_e\n",
    "    else:\n",
    "        F_e = A_e * Phi_e.T @ Phi_e\n",
    "        \n",
    "        F[np.ix_(e,e)] = F[np.ix_(e,e)] + F_e\n",
    "\n",
    "print('size: {} x {}'.format(*F.shape))\n",
    "print('rank: {}'.format(np.linalg.matrix_rank(F)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a01a3f7-c7e0-40b2-aeef-6ac3d2787a52",
   "metadata": {},
   "source": [
    "## 4. Assemble kinematic vector (with source, Dirichlet problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c855dbeb-fe00-4b49-8f23-8514e2fd3e53",
   "metadata": {},
   "outputs": [],
   "source": [
    "du = np.zeros(num_dof)\n",
    "\n",
    "##########################\n",
    "#   Dirichlet boundary   #\n",
    "##########################\n",
    "for edge in dir_edge_1:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir1, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_2:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir2, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_3:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir3, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L\n",
    "\n",
    "for edge in dir_edge_4:\n",
    "    i, j = edge\n",
    "    L = np.linalg.norm(node[i]-node[j])\n",
    "    up = gauss_line_int(dir4, *node[i], *node[j])\n",
    "    du[i] = du[i]-up/L\n",
    "    du[j] = du[j]+up/L    \n",
    "    \n",
    "##########################################\n",
    "#    Neumann boundary (without source)   #\n",
    "##########################################\n",
    "s_N = np.zeros(num_dof)\n",
    "A = np.eye(num_dof)\n",
    "\n",
    "num_red = len(ccw_hole_neu_edges)\n",
    "num_eff = num_dof - num_red\n",
    "\n",
    "# 1. Externel edge\n",
    "# effective variable (s1) is at first node of first Neumann edge \n",
    "# (edge's nodes are ordered counterclockwise within the local element).\n",
    "# ...\n",
    "\n",
    "# 2. Hole edge\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "hole_eff_dofs = [hole_neu_edge_1[0][0]]\n",
    "hole_red_dofs = [edge[1] for edge in ccw_hole_neu_edges[:-1]]\n",
    "\n",
    "A[np.ix_(hole_red_dofs, hole_eff_dofs)] = 1\n",
    "\n",
    "A = np.delete(A, hole_red_dofs, axis=1)\n",
    "\n",
    "# Neumann terms\n",
    "####### TODO : crossing edge #########\n",
    "# externel (counter-clockwise)\n",
    "# ...\n",
    "\n",
    "# hole (clockwise)\n",
    "# last edge should be excluded for closed Neumann boundary. \n",
    "fl_bar = 0\n",
    "for edge in ccw_hole_neu_edge_4:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu4, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_3:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu3, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_2:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu2, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "for edge in ccw_hole_neu_edge_1[:-1]:\n",
    "    i, j = edge\n",
    "    fl_bar = fl_bar + gauss_line_int(hole_neu1, *node[i], *node[j])\n",
    "    s_N[j] = fl_bar\n",
    "    \n",
    "\n",
    "# Source \n",
    "# Particular solution q_0\n",
    "def qx_0(x, y):\n",
    "    return x\n",
    "\n",
    "def qy_0(x, y):\n",
    "    return np.zeros_like(x)\n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    \n",
    "    Phi_e = (1/(2*A_e))*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                  [y2-y3,y3-y1,y1,-y2]])\n",
    "    \n",
    "    int_qx = gauss_tri_int(qx_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "    int_qy = gauss_tri_int(qy_0, x1, y1, x2, y2, x3, y3, degree=2)\n",
    "    \n",
    "    par = Phi_e.T @ np.array([int_qx, int_qy])\n",
    "    \n",
    "    du[e] = du[e] - par\n",
    "\n",
    "# Neumann boundary with source.\n",
    "# ..."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5619f835-d960-4905-a43c-5b77ff0dcc40",
   "metadata": {},
   "source": [
    "## 5. Constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0dba5f60-fa70-4d41-ac50-72d386d3aa69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1\n"
     ]
    }
   ],
   "source": [
    "FU = F[:]\n",
    "dU = du[:]\n",
    "\n",
    "F = A.T @ FU @ A\n",
    "d = A.T @ (dU - FU@s_N)\n",
    "\n",
    "r1 = 1\n",
    "\n",
    "print(r1)\n",
    "\n",
    "F[-1,:] = 0\n",
    "F[-1,-1] = 1\n",
    "d[-1] = r1\n",
    "\n",
    "F[0,:] = 0\n",
    "F[0,0] = 1\n",
    "d[0] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4ffc26-a9f0-408c-80eb-723329c3b8ee",
   "metadata": {},
   "source": [
    "## 6. Solve and show"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "975aa11a-a4b7-49fb-869b-39843c1242e6",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "1.0\n",
      "1.2101759607195545e-13\n"
     ]
    }
   ],
   "source": [
    "s_eff = np.linalg.inv(F) @ d\n",
    "s = A @ s_eff + s_N\n",
    "print(s_eff[-1])\n",
    "print(d[-1])\n",
    "print( np.linalg.norm(F @ s_eff - d))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c8e4ab57-f7ab-4b05-a27f-61405d7d3079",
   "metadata": {},
   "source": [
    "## 7. Construct the flux field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d4a4a3c5-081a-43ff-9f2a-1603428996bf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<visvis.wibjects.colorWibjects.Colorbar object at 0x000001E497DC41F0>\n",
      "ERROR calling '_OnAxesPositionChange':\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\wibjects\\colorWibjects.py\", line 622, in _OnAxesPositionChange\n",
      "    self.position.x = axes.position.width+5\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\misc.py\", line 217, in fsetWithDraw\n",
      "    fset(self, *args)\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 1126, in fset\n",
      "    self._Update()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 928, in _Update\n",
      "    self._CalculateInPixels()\n",
      "  File \"D:\\20220326(EFEM-Poisson)\\EFEM\\EFEM\\lib\\site-packages\\visvis\\core\\base.py\", line 997, in _CalculateInPixels\n",
      "    raise Exception(\"Can only calculate the position in pixels\"+\n",
      "Exception: Can only calculate the position in pixels if the owner has a parent!\n",
      "\n"
     ]
    }
   ],
   "source": [
    "q_h = np.zeros((E,2)) \n",
    "\n",
    "for i, e in enumerate(element):\n",
    "    x1,y1 = node[e][0]\n",
    "    x2,y2 = node[e][1]\n",
    "    x3,y3 = node[e][2]\n",
    "    \n",
    "    A_e = 1/2*np.linalg.det(np.c_[node[e], np.ones((3,1))])\n",
    "    Phi_e = 1/(2*A_e)*np.array([[x2-x3,x3-x1,x1-x2],\n",
    "                                [y2-y3,y3-y1,y1-y2]])\n",
    "    \n",
    "    if i in Eh:\n",
    "        L1 = np.sqrt((x2-x1)**2+(y2-y1)**2)\n",
    "        L2 = np.sqrt((x3-x2)**2+(y3-y2)**2)\n",
    "\n",
    "        left = np.linalg.inv([[(y2-y1)/L1,(x1-x2)/L1],\n",
    "                              [(y3-y2)/L2,(x2-x3)/L2]])\n",
    "        left = np.linalg.inv([[(y1-y2)/L1,(x2-x1)/L1],\n",
    "                              [(y2-y3)/L2,(x3-x2)/L2]])\n",
    "\n",
    "        Phi_L = np.c_[left, np.zeros((2,1))]\n",
    "        \n",
    "        Phi_e = np.c_[Phi_e, Phi_L @ Eh[i]]\n",
    "        \n",
    "        ey = np.append(e, -1)\n",
    "\n",
    "        q_e = Phi_e @ s[ey]\n",
    "    else:\n",
    "        q_e = Phi_e @ s[e]\n",
    "    # q_e = Phi_e @ s[e]    \n",
    "    q_h[i] = q_e\n",
    "\n",
    "cfv.figure()\n",
    "\n",
    "# Draw the mesh.\n",
    "\n",
    "cfv.drawElementValues(\n",
    "    ev=q_h[:,0],\n",
    "    coords=coords,\n",
    "    edof=edof,\n",
    "    dofs_per_node=mesh.dofsPerNode,\n",
    "    el_type=mesh.elType,\n",
    "    title=\"Example 01\"\n",
    "        )\n",
    "cfv.showAndWait()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "d7b83daa-f7c6-4f98-aeb6-0256c960d087",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.03504291382534708\n"
     ]
    }
   ],
   "source": [
    "fp1 = lambda x: -np.cosh(0.4)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: -np.cos(0.6)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: np.cosh(0.6)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: np.cos(0.4)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(0.6)-fp1(0.4))+(fp2(0.6)-fp2(0.4))+(fp3(0.4)-fp3(0.6))+(fp4(0.4)-fp4(0.6)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "d717b115-2bdb-4587-a873-e504a0a5460d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.1102230246251565e-16\n"
     ]
    }
   ],
   "source": [
    "fp1 = lambda x: np.cosh(0)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp2 = lambda y: np.cos(1)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "fp3 = lambda x: -np.cosh(1)/np.sin(1)/np.sinh(1)*np.cos(x)\n",
    "fp4 = lambda y: -np.cos(0)/np.sin(1)/np.sinh(1)*np.cosh(y)\n",
    "\n",
    "\n",
    "print((fp1(1)-fp1(0))+(fp2(1)-fp2(0))+(fp3(1)-fp3(0))+(fp4(1)-fp4(0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c8d29a66-eea7-4e82-b2d5-a74a26af9427",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EFEM",
   "language": "python",
   "name": "efem"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
